#include "include/Region.hpp"
#include "stdlib.h"

extern int element_shape;
extern int degree_of_basis_function;

Region::Region(std::vector<Node> n)
{     
    for(auto iter=n.begin(); iter<n.end(); iter++)
    {
        nodes.push_back(*iter);    
    }

    if(element_shape == 3)
    {
        if( degree_of_basis_function == 1)
        {
            for(int i=0; i<3; i++)
            { 
                auto basis = [i](double x, double y)
                {
                    switch(i)
                    {
                        case 0:
                        {   
                            return 1.0 - x - y;
                        }    
                        case 1:
                        {
                            return x;
                        }
                        case 2:
                        {
                            return y;
                        }
                        default:
                        {
                            std::cerr << "Invalid i=" << i << "to generate basis function";
                            abort();
                            return 0.0;
                        }
                    }
                };
                
                auto grad_x = [i](double x, double y)
                {
                    switch(i)
                    {
                        case 0:
                        {
                            return -1.0; 
                        }
                        case 1:
                        {
                            return 1.0;
                        }
                        case 2:
                        {
                            return 0.0;
                        }
                        default:
                        {
                            std::cerr << "Invalid i=" << i << " to generate basis gradient" << std::endl;
                            abort();
                            return 0.0;
                        }
                    }   
                };

                auto grad_y = [i](double x, double y)
                {
                    switch(i)
                    {
                        case 0:
                        {
                            return -1.0;
                        }
                        case 1:
                        {
                            return 0.0;
                        }
                        case 2:
                        {
                            return 1.0;
                        }
                        default:
                        {
                            std::cerr << "Invalid i=" << i << " to generate basis gradiend" << std::endl;
                            abort();
                            return 0.0;
                        }
                    }
                };
                
                basis_function.push_back(basis);
                std::vector<std::function<double(double, double)> > grad_basis;
                grad_basis.push_back(grad_x);
                grad_basis.push_back(grad_y);
                basis_gradient.push_back(grad_basis);
            }
            
        }

        if(degree_of_basis_function == 2)
        {
            for(int i=0; i<6; ++i)
            {    
                auto basis = [i](double x, double y)
                {
                    switch(i)
                    {
                        case 0:
                        {
	                        return (1.0 - x - y) * (1.0 - 2.0 * x - 2.0 * y);
                        }
                        case 1:
                        {
                            return x * (2.0 * x - 1.0);
                        }
                        case 2:
                        {
	                        return y * (2.0 * y - 1.0);
                        }
                        case 3:
                        {
                            return 4.0 * x * y;
                        }
                        case 4:
                        {
                            return 4.0 * (1.0 - x - y) * y;
                        }
                        case 5:
                        {
	                        return 4.0 * (1.0 - x - y) * x;
                        }
                    }
                };
                
                auto grad_x = [i](double x, double y)
                {
                    switch(i)
                    {
                        case 0:
                        {
                            return 4.0 * x + 4.0 * y - 3.0; 
                        }
                        case 1:
                        {
                            return 4.0 * x - 1.0;
                        }
                        case 2:
                        {
                            return 0.0;
                        }
                        case 3:
                        {
                            return 4.0 * y;
                        }
                        case 4:
                        {
                            return -4.0 * y;
                        }
                        case 5:
                        {
                            return 4.0 - 8.0 * x - 4.0 * y;
                        }
                        default:
                        {
                            std::cerr << "Invalid i="  << i << " to generate basis_gradient" << std::endl;
                            abort();
                            return 0.0;
                        }
                    }
                };

                auto grad_y = [i](double x, double y)
                {
                    switch(i)
                    {
                        case 0:
                        {
                            return 4.0 * x + 4.0 * y - 3.0; 
                        }
                        case 1:
                        {
                            return 0.0;
                        }
                        case 2:
                        {
                            return 4.0 * y - 1.0;
                        }
                        case 3:
                        {
                            return 4.0 * x;
                        }
                        case 4:
                        {
                            return 4.0 - 4.0 * x - 8.0 * y;
                        }
                        case 5:
                        {
                            return -4.0 * x;
                        }
                        default:
                        {
                            std::cerr << "Invalid i="  << i << " to generate basis_gradient" << std::endl;
                            abort();
                            return 0.0;
                        }
                    }
                };
                
                basis_function.push_back(basis);
                std::vector<std::function<double(double, double)> > grad_basis;
                grad_basis.push_back(grad_x);
                grad_basis.push_back(grad_y);
                basis_gradient.push_back(grad_basis);
            } 
        }
    }
    
    else
    {
        if(element_shape == 3 and degree_of_basis_function == 2)
        {
            std::cerr << "Sorry! Triangle elements with quadratic basis function isn't ready yet!" << std::endl;
            abort();
        }
        if(element_shape == 4)
        {
            std::cerr << "Sorry! Rectangle elements isn't ready yet!" << std::endl;
            abort();
        }
    }
    std::cerr << "Region for element shape=" << element_shape << " with " << basis_function.size() << " basis functions is ready" << std::endl;
}

std::ostream& operator<<(std::ostream& os, Region& r)
{
    os << "Reference region with following nodes' coordinates and weight " << std::endl;
    for(auto iter=r.nodes.begin(); iter<r.nodes.end(); iter++)
    {
        os << (*iter);
    }
    return os;
}

double Region::integral(std::function<double(double, double)> f)
{
    double rst = 0.0;
    for(auto iter=nodes.begin(); iter<nodes.end(); iter++)
    {
        rst += f((*iter)[0], (*iter)[1])* (*iter).get_weight();
    }
    return rst/2.0;
}
